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Abstract 


Research supported by NASA Langley Research Center includes many applications 
of aerospace design optimization and is conducted by teams of applied mathematicians and 
aerospace engineers. This paper investigates the benefits from this combined expertise in 
formulating and solving integer and combinatorial optimization problems. Applications 
range from the design of large space antennas to interior noise control. A typical problem, 
for example, seeks the optimal locations for vibration-damping devices on an orbiting 
platform and is expressed as a mixed/integer linear programming problem with more than 
1500 design variables. 

Introduction 

The purpose of this effort is to investigate the interchange of ideas between 
aerospace engineers and applied mathematicians in formulating and solving design 
optimization problems. This research also describes and provides examples of integer and 
combinatorial optimization applications that have been studied at NASA Langley Research 
Center. 

Several optimization methods, including simulated annealing, tabu search and 
branch and bound, for solving combinatorial optimization problems are mentioned herein. 
These are heuristic (i.e., rule-based) methods and are not guaranteed to converge to a 
global minimum. A brief characterization of each method is given below. For a detailed 
description, the interested reader is referred to the appropriate references. 

Metaheuristics such as simulated annealing and tabu search provide a shell within 
which a variety of other heuristics may be implemented. The definitions and notations that 
follow are taken from references 1-3. 

Let Z denote the set of all feasible states (i.e., the set of feasible solutions to the 
minimization problem generated from all possible combinations of the design variables) and 
let S denote an element of Z. To differentiate between states we define a criterion function 
c and refer to c(S) as the cost of state S. To move from one state to another we define a 
move set A and a move 8 e A. The outcome of applying all legal moves 8 e A to the 

current state S defines the set of states reachable from S; this set is typically called the 
neighborhood of S. The value of a move is the difference between the cost of the new state 

and the cost of S (i.e., c[8(S)] - c(S)). 
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Each metaheuristic begins with an initial state S 0 , chosen either at random or 
constructed algorithmically. Then the metaheuristic generates a sequence of moves {8 0 , 8j, 

8„} which determines a sequence of states through which the search proceeds. 

The mechanism by which a move is selected is one of the crucial differences 
between simulated annealing and tabu search. To appreciate the difference consider an 
improvement scheme that at state S t selects the greatest available one-move improvement. 

That is, the next move 8, is chosen based on minimum cost: 

c[8t(S t )] = nunc[8(S t )] (1) 

Equation (1) is called a greedy local improvement scheme and a metaheuristic based on 
equation (1) is called a greedy search algorithm . The greedy search terminates when no 
improving move is found. The final state of a greedy search is a local minimum (with 
respect to a particular move structure) and is, generally, not the global minimum. Both 
simulated annealing and tabu search are attempts to circumvent this difficulty. In simulated 
annealing nonimproving moves are initially accepted with a high probability; the probability 
is gradually decreased. Simple versions of tabu search attempt to avoid entrapment in local 
minima by maintaining a list of previously selected moves and deleting them from the move 

set A for a state S to avoid a return to a previously observed state. More sophisticated 

features of tabu search involve explicit use of the search history to diversify or intensify the 
search. 

The branch and bound algorithm with linear programming (LP) relaxations (e.g., 
see ref. 4) is an alternative to the above heuristic algorithms. This technique is a good 
choice for combinatorial optimization problems that involve binary design variables and 
linear criterion functions and linear constraints. The method solves a sequence of LP 
problems that establish upper and lower bounds on the solution to the integer linear 
programming (ILP) problem. These bounds are used to "prune" branches from the binary 

tree which describes the state space E. The method terminates after a fixed number of LP 

problems have been considered or when the difference between the newest upper and lower 
bounds is small compared with the modeling or measurement uncertainty. 

Three aerospace integer or combinatorial optimization problems are cited in this 
paper. The first involves the selection of the best assembly sequence for a large space 
antenna; the second seeks optimal locations for vibration-damping devices on an orbiting 
space platform; and the third involves determination of the number and position of active 
structural acoustic control actuators. For each case, the mathematical problem statement is 
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given, a solution method is suggested, and typical results are examined. Moreover, the 
contributions of mathematicians and engineers to the solution of the problem are 
acknowledged; and obstacles encountered in the solution process are reviewed. 

Antenna Assembly Sequence Problem 

Assume that an antenna (fig. 1(a)) is to be designed and erected in space using a 
large number n of truss elements. For research purposes, the antenna structure is designed 
as a tetrahedral truss (fig. 1(b)) with a flat top surface (i.e., all nodes in the top surface of 
the finite-element model are coplanar). The lengths of the truss elements must be identical 
to minimize surface distortion and to avoid the buildup of internal forces during the 
assembly process. However, because of unavoidable manufacturing limitations, the 
lengths are never precisely identical. Each truss element j has a small but measurable error 
ej . One way to minimize the impact of these errors is to assemble the antenna in such a 
way that the errors offset one another. 

A combinatorial optimization problem that determines the best arrangement of the 
truss elements is developed here. First, the objective is stated as the minimization of the 
squared L 2 norm of the surface distortion : 


d 2 = e T U T DUe = e T He (2) 

where e is the vector of measured errors; U is the influence matrix such that u tj gives the 
influence of a unit error in member j on the surface distortion at node i; and D is a positive 
semidefinite weighting matrix that denotes the relative importance of each node i at which 
distortion is measured. The matrix U can be calculated by any structural analysis software 
package, and the matrix D is often the identity matrix. With these definitions, the antenna 
distortion problem is stated as 


minimize ( 3 ) 

i~ i >=i 

over all permutations of the error vector e. Clearly equation (3) is a quadratic assignment 
problem; however this equation is not the typical version studied in the operations research 
literature. 

For the antenna assembly sequence problem, engineering insight led directly to a 
mathematical solution process. The authors of reference 5 observe that the cost of truss 
elements increases dramatically when unusual precision in length is demanded. They 
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suggest that truss elements with standard precision would be adequate if assembled in the 
correct order. Because construction of the antenna in space requires careful planning of the 
assembly sequence, the mathematical assignment of truss elements to specific antenna 
locations does not increase the complexity of the process. 

In reference 5 the following conceptual solution is proposed. The antenna is 
assembled with a random assignment of elements. Pairs of elements are selected at random 
and interchanged. If the surface distortion degrades then the interchange is reversed, 
otherwise it remains. This process continues until no further improvement is realized. The 
effect of interchanges can be predicted using equation (2) so that no hardware changes must 
be made until the best arrangement has been identified (ref. 5). 

The pairwise interchange heuristic that is suggested in reference 5 is similar to the 
greedy heuristic discussed in the introduction. This interchange heuristic is inferior to the 
simulated annealing or tabu search algorithms developed in references 6-7. For example. 
Figure 2 shows a comparison of the simulated annealing results with those of the pairwise 
interchange. In each trial both methods are initialized with the same random assembly 
sequence and the results of ten trials are plotted. Both methods result in final assembly 
sequences which are orders of magnitude better than the initial random sequence. 

However, the simulated annealing method consistently results in an excellent solution, and 
the pairwise interchange usually converges to an inferior solution. 

Suggestions from aerospace engineers for the antenna assembly sequence problem 
led directly to a good problem formulation and to a workable combinatorial optimization 
scheme. Applied mathematicians suggested an improved optimization scheme. The next 
two case studies suggest that formulating the correct optimization problem generally 
requires more interaction between mathematicians and engineers. 

Damper Placement Problem (DPP) 

One aspect of the NASA Controls-Structures Interaction (CSI) project was a set of 
laboratory experiments investigating the control of space structures. The Phase 1 CSI 
Evolutionary Model (CEM) was a large, flexible structure assembled from truss elements 
and antenna support members. (See fig. 3.) The CEM was designed to simulate some 
characteristics of a large earth-observation platform. The CEM was suspended by cables 
and was dynamically tested in the NASA Langley Space Structures laboratory. After the 
dynamic characteristics of the original model were measured, numerous active control 
concepts were applied and tested. 
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In one active control concept (ref. 8), 8 of the 1507 truss elements were removed 
and replaced by active struts. An active strut is a combination of an actuator and a sensor. 
Active struts can sense axial compression or tension and use a feedback control law to 
dissipate strain. Active struts that are placed in high-strain locations can enhance vibration 
damping. 

The goal of the DPP is to determine optimal locations for 8 active struts so as to 
maximize the minimum modal damping ratio over the first 26 characteristic vibration modes 
of the structure. As explained in reference 8, the goal is to improve damping in all target 
modes so that any vibration induced in the structure will decay quickly. In reference 8, the 
DPP is expressed as a mixed ILP problem with 1508 design variables and 27 constraints: 
maximize (5 

such that ^ v ij x j - w iP for i = 1,2,. ..26 

' ( 4 ) 

and 2^, x j -8 

j 

design variables : (5, , ..., x 15a7 

where is the fraction of axial strain energy in mode i and truss element j; p is a real- 
valued design variable that is proportional to the minimum modal damping ratio; and x is a 
vector of binary design variables such that x ; = 1 if truss element j is to be replaced. The 
optional vector of weighting factors w can be used if the control of certain modes is 
particularly important. 

For this second case study, engineering insight led to a useful description of the 
physical problem but did not provide an effective mathematical solution method. One 

possible solution method is to select the location j with the maximum value of v (j for each 

mode i. This method is reasonable for controlling one or two modes. However, if 8 struts 
must provide damping for 26 modes, then locations that simultaneously provide damping 
in several different modes must be sought. 

Reference 9 contains information used to describe and formulate the DPP. The 

authors suggest the use of the sum of axial strain energy ( v i:/ jc y ) as a measure of 
damping in mode i. In addition, they suggest that the weights should be proportional to the 

percent of modal damping in mode i (i.e., w, °c Vy ). This insight is important because 

j 

many of the first 26 modes in the CEM involve motion of suspension cables and 
deformation of the antenna support elements. The sum of modal strain energy due to 
tension or compression of the truss elements is tiny for such modes. Modes that cause little 
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or no strain in the truss elements cannot be controlled by placing active struts in truss 
locations. This phenomenon suggests that w ; .= 0 if the modal strain energy is small (e.g., 
less than 30 percent) and w~l otherwise. 

Reference 9 provides relevant information in regard to the physics of the DPP but 
does not provide efficient solution techniques. Instead of stating the problem as in equation 
(4), the authors pose the following unconstrained nonlinear programming problem: 


max 



->26 

^1507 

f X j Vi j/ 1 

1/4' 

26 1 

U=1 

_2,j=l 

l AiJJ 



(5) 


design variables : xj 


where £ is a target value of damping in mode i. The solution method is simulated 
annealing. 

Equation (4) is preferred over equation (5) for three practical reasons. First, 
equation (5) will not necessarily provide damping in each controllable mode. Second, no 
straightforward method exists for adding topological constraints to equation (5). Third, the 
effort required to solve equation (5) by simulated annealing increases with the size of the 
search domain: 

AM 

size =t — -r (6) 

[M!(/V - M)!] 


where N = 1507 is the number of possible locations and M = 8 is the number of active 
struts. The number of combinations of 1507 locations taken 8 at a time is approximately 
10 20 ; however the size increases dramatically if M increases. 

Thus, the DPP is a case in which mathematical expertise is beneficial to formulating 
the design optimization problem. For example, reference 8 discusses the solution of 
equation (4) as a mixed ILP problem. The branch and bound algorithm using linear 
programming relaxations is demonstrated. Topological constraints (e.g., a restriction on 
the selection of adjacent locations) are quite easy to add. Furthermore, the efficiency of the 
branch and bound algorithm is sensitive to the number of modes and the number of 

possible locations (i.e., the dimension of the v matrix) but not to the number of active 
struts. 

The solution to the DPP using branch and bound algorithm surprised both the 
engineers and the mathematicians. Figure 4 illustrates one solution to equation (4) for 
selecting locations of eight active struts. This solution was surprising from an engineering 
standpoint because two locations were selected on the suspension arms near the place 
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where the cables were attached to the structure. Because these arms were designed to be 
rigid supports for the flexible truss structure, they were considered unlikely locations for 

active struts. However, further analysis revealed that the v values for these locations were 

quite large in several modes. These large values of predicted axial strain energy were 
partially attributable to a modeling error and led to an improved finite element model for the 
CEM and a new set of optimal locations. However, the large values were also due to the 
basic design of the CEM. The next version of the CEM (i.e.. Phase 2) had more rigid 
support arms. Thus, engineering insights gained from solving the mathematical DPP were 
not only instrumental in finding the best locations for active struts on the Phase 1 CEM but 
also influenced the design of the Phase 2 CEM. 

On the other hand, the experimental results provided important insight to the 
mathematicians. When active struts were tested in the predicted “optimal” locations they 
provided little vibration damping for several modes. It was concluded that the method of 
finding locations using equation (4) has two weaknesses. First, the assumption is made 
that the structural finite element model is a perfect representation for the CEM. A second 
assumption asserts that the active struts are identical (i.e., they have the same mass and 
stiffness properties) to the truss elements they replace. Neither assumption is justified. 

An improved version of the DPP would include some uncertainty in specifying the 
structural finite element model. However, this uncertainty creates mathematical difficulties. 
For example, if the active struts are significantly different from the truss elements that they 
replace, then each change in the solution vector x requires a new structural model, a new 

set of characteristic modes, and a new set of v values. If the number of modes and the 

values of v are functions of the design variables x then the branch and bound solution to 

equation (4) becomes impossible and a simulated annealing approach is more appropriate. 

The DPP illustrates the need for engineering and mathematical input and the mutual 
benefits that can be gained in the optimization of engineering systems. However, important 
questions are raised in regard to the effect of both modeling errors and uncertainty on the 
optimization process. 

Active Structural Acoustic Control (ASAC) 

Assume that an aircraft fuselage is represented as a cylinder with rigid end caps 
(fig. 5) and that a propeller is represented as a point monopole with a frequency equal to 
some multiple of the blade passage frequency. Piezoelectric (PZT) actuators bonded to the 
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fuselage skin are represented as line force distributions in the x and 0 directions. Using 

this simplified model, the point monopole produces predictable pressure waves that are 
exterior to the cylinder. These periodic pressure changes cause predictable structural 
vibrations in the cylinder wall and predictable noise levels in the interior space. The interior 
noise level at any discrete microphone location can be dramatically reduced by using the 
PZT actuators to modify the vibration of the cylinder. For a given set of microphones and 
a given set of actuator locations, the control forces that minimize the L 2 norm of the noise 
are known. However, methods for choosing the optimal locations for the microphones and 
the optimal locations for the actuators have not been considered. 

The use of active structural acoustic control in cylindrical fuselage structures is 
explained in reference 10 and verified by numerous experiments (e.g., refs. 10-12). The 
results in reference 10 demonstrate that the amount of noise control depends both on the 
geometry of the source plus the cylinder system and on the locations of discrete control and 
measuring points. The force limitations of the PZT actuators must be considered in 
planning the control strategy. In addition, effective noise control strategies can either 
reduce the vibration of the cylinder or can increase the vibration of the cylinder, which 
shifts the energy to shell modes that do not couple efficiently with acoustic modes. This 
insight is important because aircraft manufacturers may reject a noise control method that 
increases vibration and in turn increases fatigue of the airframe. 

In accordance with the notation given in reference 1 1, the ASAC optimization 
problem is to minimize the sum of squared pressures at a discrete set of interior 
microphones: 

E = f. A.A'. (7) 

m= 1 

where Np is the number of microphones and * indicates the complex conjugate. The 
response at microphone m is given as: 

Nc 

K='LH mk c k +p m ( 8 ) 

*=1 

where p m is the response with no active control and is a complex- valued transfer matrix 
that represents the pressure at microphone m due to a unit control force (Iql = 1) at the PZT 
actuator k. The values in the transfer matrix can be collected experimentally (ref. 1 1) or 
they can be simulated (ref. 10). 

The cost function can be written as in equation (7) or expressed on a decibel scale 
which compares the interior pressure norm with and without ASAC: 
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For a fixed set of Nc actuators, the forces c k which minimize either equation (7) or 
equation (9) can be determined by solving a complex least-squares problem (ref. 10). 
Unfortunately, the solution vector may contain values of c k that exceed the maximum 
allowable control force. Also, the solution vector may decrease the interior noise level and 
increase the shell vibration level. (Note that an equation similar to equation (9) exists that 
compares the vibration norm with and without ASAC. A positive vibration level signifies 
an increase in shell vibration due to the action of PZT actuators.) 


In the ASAC case, engineering input complicated the optimization process. The 
engineering approach assumed that the forces c k were variable but that the actuator 
locations were fixed. Several attempts were made to use multiobjective optimization to 
trade off noise reduction and vibration reduction while imposing force constraints. These 
attempts met with limited success (ref. 1 1). The weakness of this method is that it is a 
multiobjective formulation and, thus is highly sensitive to the weights placed on each 
objective. 

An alternate way to pose the problem is to make the control forces dependent 
variables and choose the number and the locations of the actuators. Given a large number 
Nc of possible locations, the alternate procedure uses tabu search or simulated annealing to 
converge to the best subset of these locations. As each proposed subset is considered, the 
vector of control forces that minimizes E (eq. (7)) is calculated and the corresponding noise 
level (eq. (9)) is used to determine the value of the proposed move. 

For many numerical experiments with differing numbers of possible locations, 
subset sizes, source frequencies and sets of interior microphones, the same trends are 
observed. Namely, the subset of actuators that reduces interior noise also reduces cylinder 
vibration. Figure 6 shows typical results. In the figure, noise and vibration levels are 
plotted versus the tabu search iteration number. The 16 best locations are chosen from a set 
of 102 possible locations. Notice that the initial set of 16 actuators reduces the noise by 13 
dB but increases the cylinder vibration by 4 dB. However, after several iterations, both 
noise and vibration levels are reduced dramatically. By adjusting the number of actuators 
up or down from 16, the noise-reduction goals can be satisfied without an increase in 
vibration and without exceeding force capacity of the PZT actuators. 
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The best locations for PZT actuators are not intuitively obvious. For example, 
figure 7 shows the grid of 102 possible locations distributed in 6 rings of 17 locations. 

Each actuator location is specified by the (x, 0, r=a) position of its center. (Recall fig. 5.) 

The acoustic monopole is located at (x=U2, 0=0, r=1.2a) where L is the cylinder length 

and a is the cylinder radius. (The dimensions of the cylinder and the frequency of the 
source are chosen to simulate typical blade passage frequencies on commuter aircraft.) The 
shaded rectangles indicate the 16 best actuator locations. Figure 7(a) shows the best 
locations for controlling interior noise due to an acoustic monopole with a frequency of 200 
Hz. Figure 7(b) indicates the change in the best locations for an acoustic monopole with a 
frequency of 275 Hz. Notice the symmetric pattern in figure 7(a) which corresponds to a 
case in which the acoustic monopole excites one dominant interior cavity mode. Notice the 
greater complexity of the pattern in figure 7(b). Here, several cavity modes of similar 
importance are excited by the 275 Hz. monopole. 

The results in figures 6 and 7 are preliminary and are based on simulated transfer 
matrices. However, they indicate the importance of actuator location in active structural 
acoustic control. Experimental tests of the actuator placement procedure are planned. In 
these tests, the transfer matrix will be constructed using measured data and the optimal 
locations will be verified experimentally. 
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Concluding Remarks 


This paper details the complicated process by which engineering design 
optimization problems are formulated and solved. Occasionally, as with the antenna 
assembly-sequence optimization, an engineering description of a problem leads directly to a 
convenient solution method. More often, with engineering input alone, a multiobjective 
problem is described for which neither the important design variables nor the appropriate 
weighting of objectives are obvious. In addition, the design optimization problem is often 
simulated by a computer code that inadequately models the physical behavior of the system. 
These shortcomings lead to elegant mathematical solutions but meaningless optimization 
results. 

This paper illustrates the benefits of a synergistic relationship between engineering 
and mathematical experts. Mathematical expertise can be used to pose a design 
optimization problem in a less ambiguous manner. Often, mathematical experiments reveal 
useful trends that were previously unsuspected or uncover weaknesses and coding errors in 
the analysis codes. The reverse is also true; unexpected optimization results and 
experimental results can be used to improve mathematical models and to revise an 
optimization problem. 


NASA Langley Research Center 
Hampton, VA 23681-0001 
October 5, 1995 
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Figures 



(a) Antenna configuration 


(b) Finite element model 


Figure 1 . Conceptual design of a large space antenna. 
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Figure 2. Comparison of solutions found by simulated annealing and pairwise interchange 

heuristics. 
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Figure 3. CSI evolutionary model (CEM). 
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Figure 4. Optimal locations for eight active struts on CEM. 
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Figure 5. Schematic of cylinder and piezoelectric actuators. 



Iteration Number 

Figure 6. Typical performance of tabu search, which shows simultaneous reduction of 

noise and vibration. 
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(a) Frequency = 200 Hz. 

Figure 7. Optimal locations for 16 piezoelectric (PZT) actuators on cylinder. 
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